In [1]:
import numpy as np
import pandas as pd
from functools import partial
from IPython.display import HTML, display
from colorcet import fire
import datashader as ds
import dask.dataframe as dd
from bokeh.io import output_notebook
from bokeh.plotting import figure, output_notebook, show
from datashader.bokeh_ext import InteractiveImage
from datashader import transfer_functions as tf
from datashader.utils import export_image
from datashader.colors import colormap_select, Greys9, Hot, viridis, inferno
from bokeh.models.tiles import WMTSTileSource
import geoviews as gv
import holoviews as hv
from holoviews.operation.datashader import aggregate, datashade
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
Utah Motor Vehicle Collisions¶
In [2]:
df = pd.read_csv('Utah/Utah crashes_2010-2017.csv', parse_dates=[['CRASHDATE', 'CRASHTIME']])
df.head()
Out[2]:
In [3]:
df.shape
Out[3]:
In [4]:
df.iloc[1].RoadView
Out[4]:
In [5]:
df_extr = df.RoadView.str.extract('.*\?db=(.*)&project=(.*)&run=(.*)&lat=(.*)&lon=(.*)')
df_extr.columns = ["db", "project", "run", "lat", "lon"]
df = pd.concat([df, df_extr], axis=1)
In [6]:
df.rename(columns={'CRASHDATE_CRASHTIME': 'DATE'}, inplace=True)
df['dayofyear'] = df.DATE.dt.dayofyear
df['year'] = df.DATE.dt.year
df['month'] = df.DATE.dt.month
df['day'] = df.DATE.dt.day
df['hour'] = df.DATE.dt.hour
df['dayofweek'] = df.DATE.dt.dayofweek
df['incident'] = 1
df.set_index('DATE', inplace=True)
In [7]:
df['lon'] = df.lon.astype(np.float)
df['lat'] = df.lat.astype(np.float)
In [8]:
df.Latitude.min(), df.Latitude.max()
Out[8]:
In [9]:
df = df[df.Latitude>100]
In [10]:
df.shape
Out[10]:
In [11]:
df.to_pickle('df_utah.pkl')
In [12]:
import calmap
Calendar hotspots for the most severe crashes¶
In [13]:
calmap.calendarplot(df[df.CrashSeverity == 5]['incident'], linewidth=0, cmap='OrRd', fig_kws=dict(figsize=(18, 14)));
Calendar hotspots for all crashes¶
In [14]:
calmap.calendarplot(df['incident'], linewidth=0, cmap='OrRd', fig_kws=dict(figsize=(18, 14)));
We see a high number of crashes in Dec-Jan and low in Feb-Mar. Crashes on weekends are lower, escpecially on Sundays. Fridays seem to have the most number of crashes. However, the most severe crashes are not in winters, but instead in the 100 days of summer.
Weekly trend by the hour¶
In [15]:
day_week_labels = ['Mon', 'Tues', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
day_week_hr_labels = [day_week_labels[day] + '-' + str(hour)
for day in range(0,7)
for hour in range(0,25)]
df['day_week_hour'] = df.dayofweek*25 + df.hour
In [16]:
import matplotlib.pyplot as plt
fig = plt.figure()
fig.set_size_inches(50,6)
fig.add_subplot(111)
ax_day_week_hr = df['day_week_hour'].value_counts().sort_index().plot(kind='bar', color='k')
ax_day_week_hr.set_xticklabels(day_week_hr_labels);
Day of Week¶
In [17]:
ax = df['dayofweek'].value_counts().sort_index().plot(kind='bar', color='k')
ax.set_xticklabels(day_week_labels);
Month of Year¶
In [18]:
df['month'].value_counts().sort_index().plot(kind='bar', color='k')
Out[18]:
Hour of Day¶
In [19]:
df['hour'].value_counts().sort_index().plot(kind='bar', color='k')
Out[19]:
Yearly trend¶
In [20]:
df['incident'].resample('Y').sum().plot()
Out[20]:
Accidents are increasing year or year with a dip in 2014
In [21]:
fig = plt.figure(figsize=(15,5))
ax = fig.add_subplot(1, 1, 1)
ax.plot(df['incident'].resample('M').sum(), 'k', label='Total crashes')
ax.plot(df[df.CrashSeverity>3]['incident'].resample('M').sum(), 'r', label='Severe crashes')
Out[21]:
In [22]:
fig = plt.figure(figsize=(15,3))
ax1 = fig.add_subplot(1, 1, 1)
ax1.plot(df['incident'].resample('M').sum(), 'b-')
ax1.set_xlabel('years')
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('All Crashes', color='b')
ax1.tick_params('y', colors='b')
ax2 = ax1.twinx()
ax2.plot(df[df.CrashSeverity==5]['incident'].resample('M').sum(), 'r-')
ax2.set_ylabel('Severe Crashes', color='r')
ax2.tick_params('y', colors='r')
fig.tight_layout()
plt.show()
The number of crashes is increasing every year, and there is a strong seasonality, with most crashes in Nov-Dec and least in Feb-Apr. Feb 2011 has unusually low crashes or missing data (more likely). However, the most severe crashes occur in the summer months.
Accidents by the roadway junction feature¶
In [23]:
df.groupby('ROADWAYJUNCTFEATUREID')['ROADWAYJUNCTFEATUREID'].count().sort_values(ascending=False)
Out[23]:
Visualize all incidents¶
In [24]:
agg = ds.Canvas(750,750).points(df, 'lon','lat')
tf.set_background(tf.shade(agg, cmap=fire),"black")
Out[24]:
Serious crashes¶
In [25]:
agg = ds.Canvas(750,750).points(df[df['CrashSeverity']>3], 'lon','lat')
tf.set_background(tf.shade(agg, cmap=fire),"black")
Out[25]:
Contributing factors¶
In [26]:
factors = df.groupby('FIRSTHARMFULEVENTID')['FIRSTHARMFULEVENTID'].count().sort_values(ascending=False)
factors
Out[26]:
In [27]:
factors[:30].plot(kind='barh', figsize=(14, 10));
Where do these factors cause collisions?¶
In [28]:
for factor in factors.index[:10]:
print
factor_df = df[df.FIRSTHARMFULEVENTID == factor]
print('{} : {} incidents,'.format(factor, len(factor_df)))
agg = ds.Canvas(750,750).points(factor_df, 'lon','lat')
display(tf.set_background(tf.shade(agg, cmap=fire),"black"))
Visualize on map¶
In [29]:
df.loc[:, 'x'], df.loc[:, 'y'] = ds.utils.lnglat_to_meters(df.lon,df.lat)
In [30]:
hv.extension('bokeh')
url = 'https://server.arcgisonline.com/arcgis/rest/services/Canvas/World_Dark_Gray_Base/MapServer/tile/{z}/{y}/{x}'
tile_opts = dict(width=1000,height=600,xaxis=None,yaxis=None,bgcolor='black',show_grid=False)
map_tiles = gv.WMTS(url).opts(style=dict(alpha=0.5), plot=tile_opts)
points = hv.Points(df, ['x', 'y'])
map_data = datashade(points, x_sampling=1, y_sampling=1, cmap=fire, width=1000, height=600)
map_tiles * map_data
Out[30]:
When and where do collisions take place¶
In [31]:
df['hourcat'] = df['hour'].astype("category")
In [32]:
y_range = df.y.min(),df.y.max()
In [33]:
x_range = df.x.min(),df.x.max()
In [34]:
extent = x_range, y_range
plot_width = int(750)
plot_height = int(plot_width//1.2)
background = "black"
export = partial(export_image, export_path="export", background=background)
cm = partial(colormap_select, reverse=(background=="black"))
def base_plot(tools='pan,wheel_zoom,reset',plot_width=plot_width, plot_height=plot_height, **plot_args):
p = figure(tools=tools, plot_width=plot_width, plot_height=plot_height,
x_range=x_range, y_range=y_range, outline_line_color=None,
min_border=0, min_border_left=0, min_border_right=0,
min_border_top=0, min_border_bottom=0, **plot_args)
p.axis.visible = False
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
return p
options = dict(line_color=None, fill_color='blue', size=5)
In [35]:
from bokeh.models import WMTSTileSource
colors = ["#FF0000","#FF3F00","#FF7F00","#FFBF00","#FFFF00","#BFFF00","#7FFF00","#3FFF00",
"#00FF00","#00FF3F","#00FF7F","#00FFBF","#00FFFF","#00BFFF","#007FFF","#003FFF",
"#0000FF","#3F00FF","#7F00FF","#BF00FF","#FF00FF","#FF00BF","#FF007F","#FF003F",]
def colorized_images(x_range, y_range, w=plot_width, h=plot_height):
cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)
agg = cvs.points(df, 'x', 'y', ds.count_cat('hourcat'))
img = tf.shade(agg, color_key=colors)
return tf.dynspread(img, threshold=0.3, max_px=4)
p = base_plot(background_fill_color=background)
p.add_tile(WMTSTileSource(url='https://server.arcgisonline.com/arcgis/rest/services/Canvas/World_Dark_Gray_Base/MapServer/tile/{z}/{y}/{x}'))
export(colorized_images(*extent),"NYC_collission_times")
InteractiveImage(p, colorized_images)
Out[35]:
Here the order of colors is roughly red (midnight), yellow (4am), green (8am), cyan (noon), blue (4pm), purple (8pm), and back to red (since hours and colors are both cyclic).
Collisions on bike / pedestrian paths¶
In [36]:
df.ROADWAYJUNCTFEATUREID.unique()
Out[36]:
In [37]:
basemap = WMTSTileSource(url='https://server.arcgisonline.com/arcgis/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}')
def create_image_dfx(x_range, y_range, w=plot_width, h=plot_height):
cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)
agg = cvs.points(dfx, 'x','y')
img = tf.shade(agg.where(agg>np.percentile(agg,90)), cmap=inferno, how='eq_hist')
return tf.dynspread(img, threshold=0.3, max_px=4)
dfx = df[df['ROADWAYJUNCTFEATUREID']=='Bike/Ped Path Intersection']
p = base_plot()
p.add_tile(basemap)
export(create_image_dfx( *extent),"BikePed Path Intersection")
InteractiveImage(p, create_image_dfx)
Out[37]:
The most dangerous roads in Utah¶
In [38]:
streets = df.groupby('RouteName')['RouteName']
count = streets.count()
count.sort_values(ascending=False, inplace=True)
print(count[:20])
Data sourced from https://data-uplan.opendata.arcgis.com/